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We present a formalism for constructing quasi-equilibrium binary black hole initial data suitable 
for numerical evolution. We construct quasi-equilibrium models by imposing an approximate he- 
lical Killing symmetry appropriate for quasi-circular orbits. We use the sum of two Kerr-Schild 
metrics as our background metric, thereby improving on conformally flat backgrounds that do not 
accommodate rotating black holes and providing a horizon-penetrating lapse, convenient for imple- 
menting black hole excision. We set inner boundary conditions at an excision radius well inside 
the apparent horizon and construct these boundary conditions to incorporate the quasi-equilibrium 
condition and recover the solution for isolated black holes in the limit of large separation. We use 
our formalism both to generate initial data for binary black hole evolutions and to construct a crude 
quasi-equilibrium, inspiral sequence for binary black holes of flxed irreducible mass. 
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I. INTRODUCTION 

The coalescence and merger of binary black holes is ex- 
pected to be one of the primary sources of gravitational 
radiation to be detected by interferometric gravitational 
wave detectors (including the Laser Interferometer Grav- 
itational Wave Observatory, LIGO, and the Laser In- 
terferometer Space Antenna, LISA). The detection and 
interpretation of black hole mergers will be greatly fa- 
cilitated by theoretical predictions for the gravitational 
waveforms produced by these events. 

For large binary separations, post-Newtonian approx- 
imations can be used to model the binary inspiral and 
gravitational wave emission to excellent accuracy jl| . For 
small binary separations, when finite size and non-linear 
effects become more important, it is expected that nu- 
merical relativity simulations will provide the most accu- 
rate models and wave form predictions. 

Constructing numerical models of the binary inspiral 
typically proceeds in two steps (see, e.g. |^ for a recent 
review). In the first step, initial data are constructed 
by solving the constraint equations of general relativity. 
These initial data, which provide a snapshot of the gravi- 
tational fields at a certain instant of time, are not unique, 
and certain freely specifiable functions have to be chosen 
in accordance with the astrophysical situation at hand 
(see also the review Q). In the second step, the ini- 
tial data are evolved dynamically forward in time, which 
provides the subsequent binary evolution and with it the 
emitted gravitational wave signal. 

To date neither one of these two steps has been solved 
completely satisfactorily. A number of groups have con- 
structed initial data describing binary black holes in 
nearly circular orbit H H H Hll Ulia, El, and there 
have been several attempts at dynamical simulations of 
binary black holes 12, 13, 14, 15, 16). While this eflfort 
has made significant progress, the numerical modeling of 
binary black hole inspiral remains an unsolved theoretical 



problem. 

Building on the success of the BSSN formulation of the 
evolution equations of general relativity [13, uM ; ^^ have 
recently developed a code that can stably evolve single, 
rotating black holes for arbitrary long times M|. This 
code adopts a simple excision scheme 20] to remove the 
black hole deep interior and its singularities from the nu- 
merical grid. Implementation of such a scheme requires a 
coordinate system that smoothly extends into the black 
hole interior ("horizon-penetrating coordinates"), as for 
example Kerr-Schild coordinates (see ^J 22] ) . Our goal 
is to use this code for dynamical simulations of binary 
black hole systems in corotating coordinates. This re- 
quires initial data that describe binary black holes in 
quasi-circular orbit in horizon-penetrating coordinates. 

The first models of binary black holes in quasi-circular 
orbits IJ, I3 adopted the Bowen-York decomposition of 
the constraint equations 23] . When combined with max- 
imal slicing and conformal fiatness, the momentum con- 
straints can be solved analytically, and only the Hamil- 
tonian constraint needs to be solved numerically. For 
generalization to spinning black holes (compare J31) it 
may be desirable to abandon conformal fiatness. A re- 
cent spectral implementation by 11] shows that in the 
extreme mass ratio limit this approach does not recover 
the Schwarzschild test particle result, which underlines 
the need for alternative solutions. Furthermore, this ap- 
proach only provides the initial data for the gravitational 
fields, and a suitable coordinate system for a subsequent 
evolution has yet to be chosen. Clearly, it is desirable to 
choose a rotating coordinate system in which the binary 
appears approximately static (i.e. a coordinate system 
that is based on the existence of an approximate helical 
Killing vector). Such a coordinate system is constructed 
in 24]. However, this coordinate system is not horizon- 
penetrating, since the lapse is not strictly positive. This 
is undesirable for the dynamical evolution and singularity 
excision (but see ^ff\; compare [25J'). 



An alternative approach, 0, Q adopted a conformal 
thin-sandwich decomposition of the constraint equations 
[3, 123, 1231 instead of the Bowen-York decomposition. 
This approach seems more appeahng for the construc- 
tion of quasi-equihbrium initial data, since it allows for 
the explicit specification of the time derivatives of the 
conformally related metric and the trace of the extrinsic 
curvature (see also 24, 28, 29].) In addition, this ap- 
proach automatically provides a coordinate system that 
reflects quasi-equilibrium. In _6, 7], this decomposition 
was combined with the conformal-imaging approach of 
IjII . In addition to leading to some inconsistencies on the 
black hole throats (compare [SOj) this again yields a lapse 
that is not strictly positive. Attempts to combine the 
thin-sandwich decomposition with a puncture approach 
|5|, l31| fail because of mutually exclusive requirements 
between the different methods l25l. 



In this paper we borrow various ideas and approaches 
from previous investigators to construct initial data that 
are better suitable for evolution with our dynamical evo- 
lution evolution code. In particular, we adopt the thin- 
sandwich decomposition of the constraint equations to- 
gether with Kerr-Schild background data. In contrast to 
|28| we set the time derivative of the trace of the extrin- 
sic curvature to zero, which we believe will result in data 
that are closer to quasi-equilibrium. On the excision sur- 
face we impose a boundary condition that is derived from 
requiring that the time derivative of the conformal factor 
vanish there. We impose circular orbits by setting the 
ADM mass equal to the Komar mass, which is equiva- 
lent to imposing a relativistic virial theorem |a, Q (see 



also mHaj.) 



We solve these equations numerically by finite differ- 
encing in Cartesian coordinates, which leads to results 
that are less accurate than those achieved with spectral 
methods (compare ■2^), but better suited for evolutions 
with our dynamical code, which also uses finite differ- 
encing and Cartesian coordinates. The accuracy require- 
ments for initial data are much less stringent than those 
for constructing accurate quasi-equilibrium sequences, 
which are typically determined from small differences 
between large numbers. As a by-product of our calcu- 
lations, we nevertheless present a crude inspiral, quasi- 
equilibrium binary sequence. 



The paper is organized as follows. In Section II we 
review the basic equations, boundary conditions and the 
construction of quasi-circular orbits. We present numer- 
ical results in Section III, and we discuss our findings 
in Section IV. We also include several Appendices with 
specifics of our numerical implementation. 



II. BASIC EQUATIONS 

A. The Thin-Sandwich Equations 

We begin by writing the metric in the ADM form 

ds^ = -a^dt^ + -i,j{dx' + P'dt)idx^ + 13^ dt), (1) 

where a is the lapse function, /3' is the shift vector, and 
7y is the spatial metric. Throughout this paper, Latin 
indices are spatial indices and run from 1 to 3, whereas 
Greek indices are spacetime indices and run from to 3. 
The Einstein equations can then be decomposed into 
the Hamiltonian constraint Ti and the momentum con- 
straint Mi 



n = R- K,.K'^ +K'^ = 0, 



(2) 
(3) 



M' = VjK'^ - V'K = 0, 
and the evolution equations 

^a^J - -2aK,j+V^l3j+Vjl3,. (4) 

+P'VfK,,+KaVjl3' + KjiV,(3', (5) 

Here we have assumed a vacuum spacetime {Tap — 0), 
and Vi, Rij and R = j^^ Rij are the covariant derivative, 
the Ricci tensor and scalar curvature associated with the 
spatial metric 7jj . The extrinsic curvature Kij is defined 
by equation Q. 

Most decompositions of the constraint equations start 
with a York-Lichnerowicz conformal decomposition of 
the metric 



7„- = 1P 7y 



(6) 



where 4' is the conformal factor and 7.^ the conformally 
related metric [33, yj ■ -"-* is also useful to decompose the 
extrinsic curvature Kij into its trace K and a tracefree 
part A,j, 

K,j = A, + \^,jK, (7) 

and to conformally transform Aij according to 

A''^ = V'-i"i*^' (8) 



(so that Ai.j = V~^y; see |35|,|3g.) With these defini- 
tions the Hamiltonian constraint (0) becomes 

V^V' - \il>R - ^il>^K^ + \i^-'A,jA'^ = 0, (9) 

where V^ and R are the covariant derivative and scalar 
curvature associated with 7ij, and V^ = V*Vi is the 
scalar Laplace operator. 

For a complete derivation of the conformal thin- 
sandwich decomposition we refer the reader to references 
[2a. I27I ISTJ l . Here we focus on the construction of binary 



black holes in quasi-equilibrium. In a corotating coordi- 
nate system one expects the gravitational field to depend 
on time only very weakly, and it is therefore natural to 
construct initial data for which as many functions as pos- 
sible have vanishing time derivative. Within the confor- 
mal thin-sandwich formalism both the time derivative of 
the conformally related metric and the extrinsic curva- 
ture appear as freely specifiable data, and it is therefore 
both possible and natural to set 



and 



dtK^Q 



dtiti = 0. 



?t7y 



Inserting the latter into Q we obtain 



A^j 



2a 



({ipr) 



where 



{txy^ = v'x^ + v^x' 



W'^iX' 



(10) 



(11) 



(12) 



(13) 



Equation (|12() can now be inserted into the Momentum 
constraint Q, which yields 



Al/3' 



where 



{lpy^^,H^)^^av^K, 



1 ~. 



Al/3' = Vj(L/3)^^ = V'/3' + -V^Vj/?^) + R',l3\ (14) 

Finally, condition H1U|) can be inserted into the trace of 
the evolution equation Q, which, after combining with 
the Hamiltonian constraint © becomes 

V2(aV') - a UipR + ^ip^K^ 

+ l^-Uyi*^) = V'/3'V,i^, (15) 



To summarize, the thin-sandwich formalism then pro- 
vides three equations 



AL/3^-(L/?)''-'V,ln(-^ 



a 



yV'K = (16) 



VV - \^R - Tl^^'K^ + \i^~^A,jA'^ = (17) 
V2(aV)-(aV')[i^+^V'^^' 

+ U-^A,,A'A = V'/3'V,/^, (18) 



for the three unknowns a, /3' and "0. The tracefree 
part of the extrinsic curvature is related to these un- 
knowns through equation p2|) . Before the equations can 
be solved, a background geometry 7jj and a background 
trace of the extrinsic curvature K has to be chosen. 



B. Kerr-Schild Background Data 

We base our choice for the freely specifiable data on 
a superposition of two Kerr black holes in Kerr-Schild 
coordinates |2j|23,|23- A Kerr-Schild metric is given by 



Qp-u — Vl^iv ~r 



2Hl^ly, 



(19) 



where 77^^ is the Minkowski metric, and /^ is a null-vector 
with respect to both the full metric and the Minkowski 
metric, g^'^l^li, = rj^'^l^li, — 0. From the spacetime met- 
ric (|19|l the spatial metric, the lapse and the shift can be 
identified as 



lij ^ ^ij + '^Hlilj, 


(20) 


a^il + 2HlH')-^/^, 


(21) 


, 2HIH^ 


(22) 


'^ 1 -t- 2HIH* ' 



For a black hole of mass M and angular momentum Ma 
at rest at the origin, H and l^ are given by 



H = 



Mr^ 



r4 -H (a • f )2 ' 



1^ = (IJ), 



rx — a X X + (a ■ x)a/r 



with 



f2-a2 /(-2_^2)2 



-I- [a-x)'' 



1/2 



(23) 
(24) 
(25) 

(26) 



For a non-rotating black hole with a = 0, H has a pole 
at the origin, whereas for rotating black holes H has a 
ring singularity. We therefore have to excise from the 
computational domain a region enclosing the singularity. 
In this paper we adopt a non-spinning Kerr-Schild back- 
ground to describe co-rotating black hole binaries in a 
co-rotating frame. 

We want to generate initial data for a spacetime con- 
taining two black holes with background masses Ma and 
Mb, velocities va and vb, and we will assume that the 
background metric has zero spin Ma. Such initial data 
can be constructed by adopting for the freely specifiable 
background data a superposition of two Kerr-Schild co- 
ordinate systems describing two individual black holes 
[21I m. The first black hole with label A has has a 
spatial metric 



"fA ij = 5i] + "^Ha I a i I 



Aj, 



(27) 



an extrinsic curvature Kaij, a lapse ua and a shift (3\. 
The trace of the extrinsic curvature is 



Ka = 



2Ma 



T\{l + 2MA/rAf'^ 



(1 + SMaAa). 



(28) 



The second black hole has a similar set of associated 
quantities which are labeled with the letter B. 

In Section III Al we have already adopted dt'yij = and 
dfK — 0, which leaves as the freely specifiable back- 
ground quantities the background metric 7^ and the 
trace of the extrinsic curvature K^ for which we choose 
the "superpositions" 



and 



li] — Stj + 2Ha I All a j + "^Hb Ibi Ib] (29) 



K = Ka + Kb. (30) 



C. Outer Boundary Conditions 

The requirement of corotation and the conditions of 
asymptotic flatness yield boundary conditions at spatial 
infinity 



V'lr^oo = 1, 


(31) 


/3^|,^oo = n{d^)\ 


(32) 


a\r^oo = 1, 


(33) 



where 51 is the angular velocity of the corotating frame. 
Since our computational domain does not extend to spa- 
tial infinity, we have to impose approximate boundary 
conditions at a finite separation. The asymptotic behav- 
ior of the metric in a binary black hole system is similar 
to that of any rotating system in an asymptotically flat 
spacetime, including a Kerr black hole. The asymptotic 
form of a Kerr black hole tells us the form of the leading- 
order, radial fall-off term of the shift that is important 
for determining the system's angular momentum (via a 
quadrature over the extrinsic curvature in ean.l53lbelow). 
To incorporate the angular momentum of the binary in 
the outer boundary condition of the shift vector, we con- 
sider the asymptotic shift of a single, rotating Kerr-Schild 
black hole and focus on the leading terms proportional 
to the spin. In our asymptotic binary shift we "correct" 
the shift associated with the non-spinning background 
metric with terms that have the same asymptotic fall-off 
as these spin-dependent terms. A similar argument for 
choosing the form of the asymptotic shift in our binary 
has been put forward in Sec. Ill E of p7| . 
For a single Kerr-Schild black hole we have 



r = ^ + 2M^^M^±^ + 0(.-), 



py 



P' 



r^ 
2My 

2Mz 



„, ,—2My —ax ^ , 
2M ^ + 0(r 



(34) 



4M^— -f 2Af 



(4M2 + a^)z 



0{r- 



where we have assumed rotation about the z-axis, a — 
(0, 0, a). Equation (|34|l suggests the following fall-off con- 
ditions at the edge of our computational grid for the bi- 
nary: 

V'-l ~ i, (35) 



py -py 

a-l 



y_ 

X 

z 

—A- 
1 



(36) 
(37) 
(38) 

(39) 



Here the boundary condition of the shift vector consists 
of an analytic part, ^% and a higher-order part. fH^ is 
the sum of the analytic shifts from each nonspinning 
black hole (a = 0) plus (Jl x r)*, which accounts for 
all shift terms except for terms due to the orbital ro- 
tation (frame-dragging) as identified above. The form 
of the higher-order terms on the right-hand side comes 
from consideration of the way in which the system's to- 
tal angular momentum is embedded in the asymptotic 
shift, as we argued above. The coefficients for the higher- 
order terms are determined by numerically fitting to the 
data immediately interior to the boundary; they are not 
given a priori. Note that the shift in Sec IV A of ^ ex- 
hibits the same asymptotic behavior as in eqn. (|36l) - (|38|l . 
including the higher-order, fall-off terms containing the 
angular momentum data. The boundary conditions of 
lq| and ours differ only in the analytical part: the back- 
ground shift in 6] is based on isotropic coordinates, while 
in our approach it is based on Kerr-Schild coordinates. 
Apart from the background, the form for the leading- 
order terms, including the rotational terms, is the same 
in both calculations (compare eqns I3()I38I in our paper 
with eqns 161-163 in 6]). 



D. Inner Boundary Conditions 

Since the metric is singular at the center of each hole, 
some part of the black hole interior has to be excised 
from the computational domain, which introduces the 
need for inner boundary conditions. In |28J the following 
set of inner boundary conditions was adopted 

-0 = 1 all boundaries, (40a) 

^* = (i\ sphere inside hole A, (40b) 

/3* = /3^ sphere inside hole B. (40c) 

Since |23| specified the lapse as either a = ip^aACtB or 
a = '0^(aA + aB — l)j no boundary condition for the lapse 
was required. We solve equation H15|) for the lapse, and 
therefore need an addition boundary condition. 

The set of conditions (|40|l is very simple to implement, 
but does not necessarily lead to quasi-equilibrium solu- 
tions. Assuming that the black holes are equilibrium (or 
"isolated" in the language of 40j) Cook jSO] derived an 
alternative set of boundary conditions (see also |4lLl4^ '). 
Unfortunately, the resulting equations are quite compli- 
cated and difficult to implement numerically. We have 
therefore chosen to adopt an alternative set of boundary 
conditions, which is motivated by the observation that 



for corotating, quasi-equilibrium black holes in a binary 
black hole system the time derivative 



should be small 13 • 

For the lapse and shift we set 



a = QfACKB 

P' = aA/3B- 



aspA 



(41) 



(42) 
(43) 



on the inner boundaries and note that these choices re- 
duce to the correct values in the limit of infinite binary 
separation. Imposing inner boundary conditions some- 
where inside the black hole horizon implicitly assumes 
that the solution should not depend on where exactly this 
condition is imposed. This suggest that the above condi- 
tion should hold not only on the boundary itself, but, to 
a certain approximation, also in the neighborhood of the 
boundary. With our choices of the conformal metric 

7y=V^'%-^'(7.^ + 7.'^-%) (44) 

and the trace of the extrinsic curvature 



we can then compute 

Vi/3* - aK 

= Vila aPb + aBPi) - aAaB(^A 4 

= 9.(aA/3^ + aB^X) + r^-»(aA/?B + 

-aAaBiKA + Kb) 
= aA(Vf /3^ - aBi^B) + aBi^tPl 

+l3'B[diaA + aA{T^ J^ - bT' J^)] 

+(3X[^^aB+aBiT'J^-AT'Ji)] 

= /3B[i9jaA + ctAdt In Vt/tb] 



(45) 



-Kb) 
"b/Ja) 

- uaKa) 



+Pa [di^B + aBdi In Vt/ta] 
= /^BaA^j In ^/tTtaTb + PX^^Bdi In VtTtaTb 
= /3^aan(^a). (46) 

Here we have used Vi/3' — aK = as well as T^ ji = 
diluy/^, T^ ji = di\n\/j, and a = ^~^/^ for both back- 
ground black holes A and B. With 7 = V'^^Ti equation 
(|46|l can be rewritten 

V,/3' ~aK^ I3'd, H^a) = /3'9, ln(aV/v/^). (47) 

According to (|41|) we expect Vi/3* — aK to be small, so 
that (|46|l becomes 



6P'd,^p + ^pp\p 



Jl- 



di In a) = 0, 



(48) 



which provides a Neumann condition for the confor- 
mal factor on the inner boundary. Collecting the inner 
boundary conditions, we then have 



P' = aA^B + "b/^a, 
6/3* 9, V = -^P'(f^j, + d,\na) 



(49) 



In the limit of infinite separation each black hole reduces 
to an isolated Kerr-Schild black hole, which satisfies the 
above conditions. 



E. Constructing quasi-equilibrium circular orbits 
and sequences 

Solving equations (|16|l - (|18|l subject to the boundary 
conditions %\f)\ - (|39|l and 1)49(1 yields a solution describ- 
ing two black holes at a particular separation d, mass M 
and angular momentum J. Sequences of constant irre- 
ducible mass binaries in circular orbit can be constructed 
as follows (see also the flow chart in Fig.^. 

Focusing on equal-mass binaries, we first choose a value 



of the irreducible mass Mi, 



ll, which remains 



constant during the slow, adiabatic binary inspiral (see 
also 0, |53)- The irreducible mass is determined from 
the area of the black hole event horizon, but in practice 
we approximate this value by computing the area of the 
apparent horizon 



Mi, 



A 

16^ 



1/2 



(50) 



We next choose a separation d, and begin the iteration 
with a trial value of the background masses Ma = Mb = 
M, which enters the background geometry 7^ and K. 
We also choose a trial value of the orbital angular veloc- 
ity rj, which enters the orbital shift in /3* in the outer 
boundary conditions (p?^ - (ISHJ- Solving equations itT?)!) 
- 1)18(1 for these values will provide a binary that is not 
necessarily in circular orbit and does not necessarily have 
the required irreducible mass. 

To impose circular orbits, we require that the system's 
ADM mass (e.g. |B]) 

Madm - t^ / 7™7^"(7mn,, - 7,«,m)d'^., (51) 
be equal to its Komar mass |43 



Mk 



1 

47r 



7*^(V,a-/?'=X,fc)d2^,. 



(52) 



a 



aA^B, 



In the above expressions d^Si — {\ /2)^^/'^ eij^dx^ dx'^ is 
the two-dimensional surface area element. In many cases, 
P'^Kik falls off faster than 0{r~^) in ((5^ so that the sec- 
ond term vanishes; in our case, however, this term cannot 
be neglected. We evaluate these integrals as described in 
Appendices El and 

The equality of the ADM and Komar masses is closely 
related to a relativistic virial theorern |38l and indicates 
that the sj^acetime is stationary p4l |45| in the rotating 
frame. In 7] this criterion was adopted to impose circular 
orbits in binary black hole spacetimes (see also |[32j for 



choose Mir 



choose separation d 



choose parameter M 



setup 7ij and K 
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choose angular velocity Q, 
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Yes 



obtain the quasi-equilibrium solution 



FIG. 1: Flow chart for the construction of sequences of quasi-equilibrium, circular-orbit binary black holes 



a pedagogical illustration). In our code we iterate over 
Q. until Madm = Afx has been achieved to within an 
accuracy of 1 part in 10^. 

For a given circular orbit we then iterate over the back- 
ground mass M until the irreducible mass Mi„ has taken 
the desired value to within an accuracy of 1 part in 10® 
|53j. Finally we vary the binary separation d to con- 
struct an approximate inspiral sequence. For each model 
the ADM mass Af adm is found from H51|l and the angular 
momentum from 



1 

8^' 



J^ = 7;^e^/ i> x^K^k(fSi 



(53) 



(see Appendix Q • The innermost stable circular orbit 
(ISCO) can be found by locating minima in the ADM 
mass and the angular momentum along a constant-mass 
sequence. 



III. NUMERICAL RESULTS 

A. Tests 

In this section we present two tests of our code. The 
first test shows second-order convergence to the analytic 
Kerr-Schild spacetime for a single black hole. The sec- 
ond test shows second-order convergence in a particular 
(non-equilibrium) binary solution previously considered 
by Pfeiffer ^. 



1. Single Black Hole Tests 

To recover a single non-rotating black hole in Kerr- 
Schild coordinates we set Mb = and located the black 
hole A at the origin. We then solve equations H16|) - H18|) 



i2^ 



i^ 
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FIG. 2: Numerical errors in the conformal factor tp and the shift /?' for three different resolutions Aa; for a single, non-rotating 
Kerr-Schild black hole. 



imposing Dirichlet inner boundary conditions 



^ 
P' 



1, 

1 
VI + 2HlHi 



(54) 



on a sphere of radius ^excision = I.SMq, where A/q is the 
background mass of one of our black holes at infinite 
separation. We run these tests with the IBM pSeries 690 
machine in NCSA. A typical run with 160 x 80 x 80 grid- 
points takes about 16 cpu hours and uses about two giga 
bytes of memory. In Fig. |21 we show numerical errors 
for the conformal factor ij) and the shift /3* for several 
different resolutions. The errors scale as expected, estab- 
lishing second-order convergence of our code. We chose 
to impose the above Dirichlet condition for ij) instead of 
the inner boundary condition H49|) because, in our nu- 



merical implementations, Dirichlet conditions are easier 
to impose at a fixed physical location, which is manda- 
tory for achieving second-order convergence. For a res- 
olution of Ax /Mo = 1/6 and with the outer boundary 
at 2GA'fo, we find J/M^ < lO'^, Madm/Mq = 0.9990, 
and Mah = 1.0013, which shows that the error in our 
solution is of the order of a fraction of a percent. 



2. Comparison with a Previous Binary Black Hole 
Calculation 

In this Section we compare with numerical results by 
Pfeiffer, Cook and Teukolsky ^ (hereafter PCT), who 
used spectral methods to construct black hole binaries. 
PCT solve equations H16|l and (|18|l . but instead of solv- 
ing H17|) for the lapse they experiment with two different 
choices for an analytic densitized lapse function. Here we 
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FIG. 3: Numerical convergence of the conformal factor ip and the shift /3* for the PCT 28| binary configuration described in the 
text. We plot the rescaled differences between four different resolutions to establish second-order convergence (see Appendix 
El. 



compare with their results for 



ai^-^ 



OtAOiB- 



(55) 



For this test we also adopt their boundary conditions 
(|40|l . The implementation of PCT allows for imposing 
the outer boundary conditions H31|) at a much larger 
separation than we can afford. Following PCT we con- 
struct a nonspinning black hole binary with background 
mass M — Afo, centers of the excised spheres at a co- 
ordinate separation d = 11.75Mo, and angular velocity 
r2 = 0.0421/Afo- (According to the effective potential 
method in jj| this choice of il corresponds to a circular 
orbit.) 

In Fig. Owe compare our results for the conformal fac- 
tor ip and the shift /3* for different resolutions (all with 
^'excision = 2.0Mo), and again establish second-order con- 
vergence or our numerical code. In Table Q] we tabulate 



both PCT's and our results for the binaries irreducible 
mass, the ADM mass, the angular momentum and the 
proper separation between the horizons. To compute the 
ADM mass and the angular momentum we extrapolate 
our numerical results to a grid with outer boundaries at 
150Mo as explained in AppendicesEl^ndini We find val- 
ues that are very similar to those found by PCT, but ours 
do not converge to theirs for a fixed location of the outer 
boundary. We believe that this is caused by the proxim- 
ity of our outer boundary, and expect, as the results in 
Table m suggest, that the agreement would improve by in- 
creasing both the resolution and the distance to the outer 
boundaries. From Tabled we estimate that the numer- 
ical errors in our results, given the grid resolution and 
location of the outer boundaries adopted in this paper, 
are of the order of about a percent. 



Aa;/M 


Domain 


Mi„/M 


Madm/M 


J/M' 


£/M 


CTS 


1.06528 


2.08436 


3.3790 


10.3971 


1/2 


small 


1.089 


2.106 


3.389 


10.121 


1/3 


small 


1.081 


2.110 


3.391 


10.123 


1/4 


small 


1.079 


2.113 


3.391 


10.125 


1/6 


small 


1.074 


2.120 


3.392 


10.126 


1/2 


big 


1.088 


2.092 


3.379 


10.122 


1/3 


big 


1.081 


2.093 


3.381 


10.124 


1/4 


big 


1.079 


2.094 


3.382 


10.126 


1/6 


big 


1.074 


2.097 


3.382 


10.127 



TABLE I: Comparison of our results with those of Pfeiffer, 
Cook and Teukolsky 12^ for a binary with M = Mq, d — 
11.75Mo and fi = 0.0421/Mo. Here M is the mass parameter 
of both black holes, Mi„ is the irreducible mass, Madm is the 
ADM mass, J is the angular momentum, and i is the proper 
distance between black hole horizons. Here a small boundary 
box refers to -20Af < x < 20M, < y < 20M and < 
z < 20M; a large boundary box refers to — 30A/ < x < 30M, 
<y < 30M and < 2 < 30M. 



Reference 


Et/fi 


J/(2^Afi„) 2nMi,, 


Schwarzschild 


-0.0572 


3.464 


0.068 


m 


-0.09030 


2.976 


0.172 


[5J 


-0.092 


2.95 


0.18 


m 


-0.068 


3.36 


0.103 


[42] 


-0.058 


3.45 


0.085 


This work 


-0.06 


3 


0.08 


[54] 


-0.0668 


3.27 


0.0883 



TABLE II: Values for the binding energy Et/n, the angular 
velocity 2nMirr and the angular momentum J / {2fiA'Iiri) at the 
ISCO as obtained in different approaches. The Schwarzschild 
results refer to the innermost stable circular orbit of a test 
particle in a Schwarzschild spacetime. Our results for this 
work are determined from the turning point of the binding 
energy curve in Fig.|3 but, as discussed in the text, they are 
prone to large errors. 



B. An approximate inspiral sequence 

We now proceed to construct an approximate inspi- 
ral sequence for a nonspinning black hole binary system, 
adopting the inner boundary conditions described in Sec- 
tion ^TOJ Contours of the conformal factor ip, the lapse 
a and the shift (3^ for one particular binary separation 
are shown in Fig. 01 

The binding energy of an equal mass binary can be 
defined as B 



Eb = Af A 



DM 



2Mu 



(56) 



A simultaneous turning point in the binding energy and 
the angular momentum locates the ISCO. In Fig. we 
show both the binding energy and the angular momen- 
tum. We show numerical results obtained with a reso- 
lution of Ax = Mq/4: and with the outer boundary at 
30Afo, which corresponds to the highest resolution and 
largest grid run in our comparison in Section flll A 21 The 



code is run with the IBM pSeries 690 machine in NCSA. 
Each run with 240 x 120 x 120 gridpoints takes about 
1600 cpu hours and uses about 8 giga bytes of memory. 
We set r-oxcision = 1.6Afo for these models. We also com- 
pare our results with the numerical results of: Cook [431 , 
for an Eddington-Finkelstein slicing and daip/dr = in- 
ner boundary condition ([43|a), an Eddington-Finkelstein 
slicing and a4> = 1/2 inner boundary condition ( 42]b), 
and a maximal slicing and datp/dr = aijj/2r inner bound- 
ary condition ('|42l|c): with the binary initial data in n; 
with the second-order, post-Newtonian sequence in jj]; 
and with the third-order, post-Newtonian sequence in 

Our results for the binding energy agree fairly well with 
those of [431 b, while our results for the angular momen- 
tum do not. Note that we find a turning point in the 
binding energy, but not in the angular momentum. There 
are several possible reasons for these findings. Solving the 
constraints in the thin-sandwich decomposition leads to 
configurations that are in approximate equilibrium, but 
lacking dynamical evolutions, it is difficult to determine 
just how good this approximation is for this scenario .5^. 
Another potential reason for our findings are the inner 
boundary conditions l|49() . which may lead to undesir- 
able deviations from quasi-equilibrium. Probably more 
important, however, is the limited numerical accuracy of 
our finite-difference, Cartesian code. From Table HI we 
estimate that the accuracy of our values for the masses 
and angular momenta is of the order of a percent or so. 
From equation H56(l the binding energy is computed as 
the difference between two masses, and is of the order 
of about 10 % of each of those masses (see Fig. O. The 
relative error in this smaller difference is therefore signif- 
icantly larger than the error in each of the terms sepa- 
rately, and may be as large as 10 % or more. Such an 
error is large enough to spoil the accuracy of an inspiral 
sequence. However, if taken at face value, the orbital pa- 
rameters at the turning point of the binding energy agree 
fairly well with recent results for the binary black hole 
ISCO (see Table |nj. 

While our results are not accurate enough to reliably 
locate the ISCO, we do believe that they are suitable 
for adoption as initial data in current dynamical evolu- 
tion calculations in finite difference implementations. For 
these purposes, the accuracy of the initial data only needs 
to be as small as that of the dynamical evolution itself. 
The individual metric quantities that must be specified 
as initial data are not small differences of large numbers 
and are determined to about a percent. 

We also note that solving the constraints in the Bowen- 
York formalism leads to higher accuracy solutions even in 
finite difference implementations (compare j^]). There, 
the momentum constraint can be solved analytically, 
leaving only the Hamiltonian constraint to be solved nu- 
merically. Moreover, the angular momentum can be de- 
termined analytically in terms of the background quan- 
tities. The Bowen-York formalism also adopts maxi- 
mal slicing, i^ = 0, so that octant symmetry can be 
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adopted and a higher grid resolution can be chosen (com- 
pare Appendix IE|l . In our approach, five coupled equa- 
tions are solved simultaneously, and all orbital parame- 
ters are computed numerically from the solutions, which 
will clearly lead to a larger numerical error. 



IV. DISCUSSION 

We present a method for constructing solutions to 
the constraint equations of general relativity, describing 
quasi-equilibrium binary black holes in nearly circular 
orbit. We expect that these solutions are suitable initial 
data for dynamical evolution with current finite differ- 
ence evolution codes. 

We solve the constraint equations in a conformal thin- 
sandwich decomposition (e.g. [23|), and impose quasi- 
circular orbits by imposing that the ADM mass of the bi- 
nary be equal to its Komar mass (compare '^). We adopt 
a superposition of two Schwarzschild black holes in Kerr- 
Schild coordinates as the conformal background solution. 
This background choice leads to horizon-penetrating co- 
ordinates, which are needed for dynamical evolutions, 
and is likely to produce less spurious gravitational radia- 
tion, at least for rotating black holes. We present a new 
set of simple inner boundary conditions, to be imposed 
on the excision surface inside the black hole, which we 
hope leads to reasonable approximation to equilibrium 
(compare j3(|). 

We present two numerical tests - one for the limiting 
case of an isolated black hole, and the other for a binary 
configuration considered in g8|. We also construct an 
approximate inspiral sequence. Our numerical accuracy 
may not be sufhcient to track the binding energy to high 
precision, since it is computed as the small difference be- 
tween two significantly larger numbers. However, we do 
expect that our solutions provide adequate initial data 
for current finite-difference evolution codes. We also ex- 
pect that when our formalism is implemented with higher 
resolution and/or more accurate numerical schemes, the 
inspiral sequence may provide a more reliable estimate of 
the ISCO parameters. 
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APPENDIX A: THE ADM MASS INTEGRAL 

Having conformally decomposed the spatial metric 7^ , 
the ADM mass integral (|51|) can be written as 



Madm 



IGtt 



(r* _ YJ'j - 8V'ip)dS, 



(Al) 



where F* = 7 F|^ (compare [19|) and where Soo is a 
closed surface at spatial infinity. For the background 
data described in Section IIIBI the first two terms can 
be evaluated analytically and yield the sum of the two 
background masses Ma + Mb, or, for equal mass bina- 
ries, 2Mo. The ADM mass therefore reduces to 



M, 



ADM 



2Mn 



1 

2^ 



W'tpdS,. 



(A2) 



We now convert this surface integral into a volume in- 
tegral using Gauss' law; volume integrals are typically 
more accurate numerically than surface integrals. How- 
ever, since a volume, say Vi, containing a black hole sin- 
gularity is excised, a surface integral over the surface of 
that volume, say Si, remains 



-Madm = 



1 



2Mo 


~ 2^ 


1 


V'lpdSi 


1 


2Mo 


1 

~ 2^ 


Jsi 


S/'ipdSi 




1 
~2^ 




^d^x 


2A/o 


1 

"2^ 


JSi 


V'^dS, 




+ 16 / ^[- 


-^R 



+iP'''A,jA'^ - -V'^i^^ld^x, 



(A3) 



Here we have used the Hamiltonian constraint in the 
third equality, and have denoted the volume outside Vi as 
V2 + ^3 + Voo as illustrated in Fig. Volume V2 denotes 
the space covered by our computational grid. Given our 
constraints on numerical grid resources, this volume ex- 
tends only to a separation of typically 30Mo from the 
black holes. Restricting the ADM integral (|A3I) to this 
volume would introduce a fairly large error. We therefore 
extend the integration to a larger volume V3 , in which the 
integrand is estimated by extrapolating (3^, a, and ip from 
their values and fall-off conditions on the outer boundary 
of the computational grid ^2 : 



^ p 


^1 + ^, 

r 


a p 


^1 + ^, 
r 


/?- p 


J.6 


py ^ 





(A4) 
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"ip contour plot 



a contour and shift vector plot 




FIG. 4: The contour plots of the binary initial data with coordinate separation 12Mo. The corresponding orbital velocity 
2Q.Mirr ~ 0.0788. The left panel shows the contour of the conformal factor tp; the right panel shows the contour of lapse a 
and the shift vector in the equatorial plane. The excision radius is I.6M0. The two thick circles in the plots are the apparent 
horizons. The numerical value of the lapse a on the apparent horizons exceeds 0.5 and is therefore positive, as required for 
horizon-penetrating coordinates. It is evident from the plots that the solutions do not satisfy octant symmetry (see Appendix 



(3^ 






Here the Ui are coefficients that are determined as follows. 
For any point in V3, say r, we find the intersection of the 
position vector r with ^2 . The value of the function at 
that intersection determines the coefficient a^. Once the 
coefficients a^ have been found, the functions ■0, ck and /3' 
and hence the integrand of the ADM mass can be evalu- 
ated in V3 (compare (3)- Typically, the boundary of V3 
is at a separation of ISOMg from the black holes, so that 
this construction increases the volume of our integration 
by a factor of 125. 



APPENDIX B: THE KOMAR MASS INTEGRAL 

The Komar mass can be defined for stationary, asymp- 
totically flat spacetimes. Stationarity implies the exis- 
tence of a Kilfing vector ^'', which can be written as 



r-an"+/3", 



(Bl) 



where n^ = — at.^t is the time-like unit normal on the 
spatial hypersurfaces Ef . The Komar mass is defined as 



87r 



47r 



(B2) 
where Soa is a closed hypersurface of Sj, diffeomorphic 
to a 2-sphere, at spatial infinity, and where we have used 



Kilhng's equation ^[^'''l = S^^'^''' . The bi-vector dS^^ = 
2nr^dS'i,i, where dSu is a spatial oriented surface area 
element, is normal on both 5oo and Ej. From IjBip we 
find 



e-'^^n. 



-a'" + /S'^'^'n. 



-a'" 



P^n^r 



(B3) 



Using the identity n^'" = —K^'^ — a^'ny , where a^ = 
n'^n^-i, is the 4-acceleration of normal observers, we ob- 
tain 



Inserting this into l|B2|) yields 



Mk 



1 

47r 



{\/'a~l3^K/)dS^ 



(B4) 



(B5) 



The term P^ K/ often falls off faster than 1/r^ in an 
asymptotically flat space, in which case its contribution 
to the integral vanishes. Here, however, this term must 
be retained. 

The Komar mass is independent of the surface S on 
which the integral is evaluated, as long as all matter 
sources are inside of S. To demonstrate this, we con- 
vert the surface integral in (jBSp into the volume integral 



Mk = ^ I {V^V'a-l3^V'K,j-K,,^'(i^)^d^x. (B6) 
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0.1 




0.1 

2^M. 



FIG. 5: The binding energy Ei, and angular momentum J as a function of orbital angular velocity. The dotted lines are least 
square fits to our numerical results, marked by open circles. We compare with thin-sandwich results in [4^ for an Eddington- 
Finkelstein slicing and datp/dr = inner boundary condition (|42|a), an Eddington-Finkelstein slicin g a nd aij) = 1/2 inner 
boundary condition ( 43lb), , and a maximal slicing and daip/dr = a'ibj2r inner boundary condition (|43]c), with the binary 
initial data in |3|, with the second-order, post-Newtonian sequence in Q], and with the third-order, post-Newtonian sequence 
in |29l . Here /x is the reduced mass (Mirr/2) and Mirr is the irreducible mass for one black hole. 




FIG. 6: The diagram illustrates the relation between the vol- 
umes Vi, V2, and V3 and the surfaces iSi, S2, and Ss,. Si is a 
boundary outside both black holes, but well inside the com- 
putational domain, 1S2 is a boundary near the outer edge of 
the computational domain, and 1S3 is a boundary well outside 
the computational domain (at very large, finite radius). 



the evolution equation Q 
V.V^a = a Ik,jK'' + ^(p + s) ) +l3'^,K-dtK, (B7) 



and the momentum constraint 



(B8) 



where, for completeness, we have included the matter 
sources p, Si and s 

s. = -l^^.n,T^''' , (B9) 

The volume integral ljB6p then becomes 

-dtK-K,jV'l3^ -I3\s,y (BIO) 

We now use the evolution equation Q) to rewrite the 
term Kij\7^f3^ as 



K.jV'P^ = aK.jK'^ + ^K'^da^j. 



(Bll) 



This brings the integral (jBlOp into the form 



This integral can be rewritten by inserting the trace of 



-dtK-^-K^'darj). 



(B12) 
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As in the calculation of the ADM mass, part of the nu- 
merical grid may have to be excluded from the integra- 
tion, for example if it contains a black hole singularity. 
The integral over an outer surface 5oo can then be writ- 
ten as a volume integral V and a surface integral over an 
inner surface, e.g. Si as in FigEl 



APPENDIX D: 



Mk = 



1 

47r 

1 

An 



{V'a - P^K/)dS^ 



V2+V3 + V, 



^d^x(_^-a{p + s)-(3's^ 



-^tK~-K'^^a^J 



An 



{V'a-P^K/)dS,. 



(B13) 



From the above assumption of stationarity, the time 
derivatives of 7^ and K have to vanish, and as long as 
there are no matter sources in fl, p = s — Si — 0, the 
volume integral vanishes and we have 



(compare [STf'). 



(B14) 



APPENDIX C: THE ANGULAR MOMENTUM 
INTEGRAL 

In Cartesian coordinates, the angular momentum can 
be defined as 



•^^ " 8^'*- 



x^K'kd^Sf, 



(CI) 



(see '23', H&i). In this paper we only consider rota- 
tions about the z-axis, and therefore compute only the 
z-component of the angular momentum which can be 
rewritten as 



J. = 



1 

1 

+ 8^ 

1 
~2 



{xA\-yA',)dSi, 



(xA'-yA'jdSi 



xA.dyf' - Ay^ 



A\ + -^'xVyK 



rryV^K 



+lyA^J^.r)V^d' 



(C2) 



As in the calculation of the ADM mass f Appendix IX|) we 
have converted the surface integral into a volume integral 
for greater numerical accuracy. As before, we evaluate 
the integral from the numerical data in volume V2 and 
from extrapolated values in volume V3. We neglect only 
those contributions to the integral from volume Vac- 



INERTIAL AND ROTATING 
FRAMES 



Rotating frames are not asymptotically flat, so that 
the expressions for the ADM (Appendix 0, angular mo- 
mentum (Appendix [nj and Komar mass f Appendix IB|) 
have to be re-evaluated. 

The barred coordinates t, x, y and z in an incrtial 
frame are related to the unbarred coordinates t, x, y and 
z in a rotating frame by the transformation 



t = t, 

X = a;cos(cji) — j/sin(cL'i), 

y = a;sin(a;t) -I- y cos(ci;i), 

z ~ z. 



(Dl) 



Here we are assuming a constant angular velocity uj — 
(0,0,0;) and rotation about the z-axis. At an arbitrary 
instant t = t = at which the two frames are aligned the 
gravitational field variables are related by j5q 



a, 



13' = P' + {uJxr)\ 









(D2) 



The only effect of this transformation is therefore the 
appearance of a new term in the shift. Since the shift 
does not enter the integrals for the ADM mass nor the 
angular momentum, those quantities remain unchanged 
and we only have to re-evaluate the Komar mass. 

Transforming between the rotating and incrtial frame, 
we find that the Komar mass in the rotating frame Mk 
is related to that in the incrtial frame Mk by 



M, 



K 



1 

An 
1 

An 



{V'a-p^K/)dS^ 
{V'a - P^K/)dS^ 



An 



UJ 



MK--—e^.tfx'K/dS, 
Mk - 2ujJ 



{lu X r)^K/dSi 



(D3) 



(note that the angular momentum is the same in both 
frames, so J = J). 



APPENDIX E: SYMMETRIES OF THE 
COMPUTATIONAL DOMAIN 

Symmetries can be used to reduce the size of the com- 
putational grid, making it desirable to incorporate as 
many symmetries as possible. One might expect that 
the binary black hole configuration studied in this pa- 
per would allow for octant symmetry. In this Appendix 
we show that this is not the case, due to the presence 
of a non-zero trace of the extrinsic curvature (i.e. non- 
maximal slicing). 
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Consider, for example, the momentum constraint Q 
and assume, for simplicity, conformal flatness. On the 
a; = plane, it would be natural to assume that /3^ be 
symmetric, while /3^ and (3^ be anti-symmetric (compare, 
for example, |33)- Computing A^^ from the shift accord- 
ing to (|12|l . assuming that all scalars are symmetric on 
all coordinate planes, shows that A^^ , A^^ , A^' and A^^ 
are all anti-symmetric on the x = plane, while A^^ and 
A'^ are symmetric. The divergence V^A^*, for example, 
is then symmetric on the x = Q plane. The gradient 
V^K, however, must be anti-symmetric, meaning that 
the momentum constraint Q violates this symmetry as- 
sumption. Similar arguments hold on the y = plane. 
Under the assumption of maximal slicing K — Q, octant 
symmetry can be adopted, but in this paper we adopt 
a Kerr-Schild background with K ^ Q. The above issue 
does not apply on the z = plane, so that equatorial 
symmetry can be assumed even in the non-maximal slic- 
ing case K ^Q. 

It is possible, however, to adopt 7r-symmetry, whereby 



times and showing that numerical errors scale in the 
expected way. Given the constraints of computational 
resources it is often impossible to double the grid size 
several times, so instead we establish second-order con- 
vergence by considering three arbitrary (but different) 
grid spacings hi, ft,2, and /13. 

Let Q{h) be a quantity obtained from a finite difference 
scheme with spacing h. A Taylor expansion around h — 
yields 



Qih) = g(o) 

= Qo + Qih 



, dQ 



h=0 
2 



2 



Q'lh 



2 dK'^ 
0[h^). 



h=0 



0(/i3) 
(Fl) 



fi-x,-y,z) = af{x,y,z). 



(El) 



Second-order convergence implies that Qi = 0. Given 
two different resolutions hi and /i2 we can eliminate Qo 
and find 



For scalar functions we have cr = 1, while for the x, y 
and z-component of the shift we have CTj; = — 1, CTy = — 1, 
and CTz = 1, respectively. 



Q2 = 



Q{h2)~Q{hi) 
hi -hi 



0{h). 



(F2) 



APPENDIX F: SECOND-ORDER 
CONVERGENCE 

Second-order convergence is most easily demonstrated 
by doubling the computational grid resolution several 



Alternatively, Q2 can be computed from the two grid 
spacings /12 and /13. Second-order convergence can then 
be established by showing that the differences between 
different values for Q2 decrease at least as fast as h. 
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